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ABSTRACT 

A new multi-dimensional Hierarchical Structure Finder (HSF) to study the phase- 
space structure of dark matter in A''-body cosmological simulations is presented. The 
algorithm depends mainly on two parameters, which control the level of connectivity of 
the detected structures and their significance compared to Poisson noise. By working 
in 6D phase-space, where contrasts are much more pronounced than in 3D position 
space, our HSF algorithm is capable of detecting subhaloes including their tidal tails, 
and can recognise other phase-space structures such as pure streams and candidate 
caustics. 

If an additional unbinding criterion is added, the algorithm can be used as a 
self-consistent halo and subhalo finder. As a test, we apply it to a large halo of the 
Millennium Simulation, where 19% of the halo mass are found to belong to bound 
substructures, which is more than what is detected with conventional 3D substructure 
finders, and an additional 23 — 36% of the total mass belongs to unbound HSF struc- 
tures. The distribution of identified phase-space density peaks is clearly bimodal: high 
peaks are dominated by the bound structures and show a small spread in their height 
distribution; low peaks belong mostly to tidal streams, as expected. However, the pro- 
jected (3D) density distribution of the structures shows that some of the streams can 
have comparable density to the bound structures in position space. 

In order to better understand what HSF provides, we examine the time evolution 
of structures, based on the merger tree history. Given the resolution limit of the 
Millennium Simulation, bound structures typically make only up to 6 orbits inside 
the main halo. The number of orbits scales approximately linearly with the redshift 
corresponding to the moment of merging of the structures with the halo. At fixed 
redshift, the larger the initial mass of the structure which enters the main halo, the 
faster it loses mass. The difference in the mass loss rate between the largest structures 
and the smallest ones can reach up to 20%. Still, HSF can identify at the present 
time at least 80% of the original content of structures with a redshift of infall as high 
as z < 0.3, which illustrates the significant power of this tool to perform dynamical 
analyses in phase-space. 

Key words: methods: data analysis, methods: numerical, galaxies: haloes, galaxies: 
structure, cosmology: dark matter 



1 INTRODUCTION 



When IZwickvl (|l933l ) studied galaxy velocities in clusters, 
he was the first to notice that there should be about 
one order of magnitude more matter in the universe than 
the observed amount of baryonic matter to explain the 
proper motions of galaxies through gravitational forces. 
Dark Matter (DM) was introduced to overcome this prob- 
lem. Later on, the existence of a dark matter component 
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was confirmed by th e analysis of galaxy rotation curves 
(iRubin fc Fordlll970D. Recent studie s of gravitational lens- 



ing (e.g. Van Waerbeke et al.l I2OO0I ') and, more generally. 



multiwavelength observations in e.g. the COSMOS project 
(Masscv ot al. 2007) provide additional proofs for the exis- 
tence of DM. Other constraints on the non-baryonic nature 
of DM were also s et by the analysis of t he Cosmic Microwave 
Background (e.g. iHinshaw et al.ll200d ). 

For the last three decades, the DM paradigm has been 
studied extensively in the context of cosmological N-body 
simulations. The comparison of structures formed in such 
simulations to observed ones excluded some of the theo- 
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retical models, such as Hot Dark Matter and led to the 
nowadays commonly accepted A-Cold Dark Matter (ACDM) 
model. In the ACDM model, dark matter is coUisionless, 
with a very small velocity dispersion at high redshift; struc- 
tures axe built in a hierarchical, bottom-up process, where 
small structures arise first, seeded from initial fluctuations, 
and then merge together to build up larger and larger struc- 
tures, designated commonly as haloes. Inside the gravita- 
tional wells of these dark matter h aloes, baryonic matter 
forms galaxies l|White fc Reeslll979l ). 

Recently, the efforts to finally identify the physical na- 
ture of dark matter particles, either directly through de- 
tecting them in ground-based dark matter particle detec- 
tors or indirectly by observing their annihilation radia- 
tion, have intensified. At the same time, th e ever increas- 
ing r esolution of N-body simulations (e.g. ISpringel et al.l 
I2OO8I ) puts new levels of demand on the field of the- 
oretical study of non-linear halos. The careful analy- 
sis of cosmological structures moves from the study 
of spherically averaged three-d imensional density profiles 
iNavarro. Frenk fc Wliit3ll997l ') to the study of the full six 
dimensional phase-space. For example, this concerns the 
analyses of the properties of caustics described analytically 
in e.g . Bertschinger's secondary infall model ( B crtschingen 
1 19851 ) and recently reviewed in the c o ntext of numerical sim- 
ulatio ns (jMohavaee fc Salatil l2008l : IWhite fc VogelsbergeJ 
I2OO8I ). Investigations of full phase-space structures in- 
clude accurat e simulati ons of two dim ensional phase-space 
ijAlard fc Colombi 20051 : IColombifcTouma, 2007, ') and anal- 
yses relying on a metric approach to six dimensional phase- 
space in A^-body sim ulations jVogelsberger et al] I2OO8I : 
IWhite fc Vogelsberged[2008 ). 



A particularly important step in understanding 
dark matter clustering lies in an analysis of the 
bound structures found in A'^-body simulations. This is 
at present usually carried out with s tructure finders 
such as SUBFIND llSpringel et"al] l200ll) . ADAPTAHOP 



ngel 

jAubert. Pichon fc Colombil bOol Tor PSB (|Kim fc ParkI 



I2OO6I ). Following this path, we present a new multi- 
dimensional Hierarchical Structure Finder which comple- 
ments all the above numerical methods with an effective and 
robust analysis of phase-space structures in full 6D space. 

The paper is organised as follows. First, we review cur- 
rent structure finders in Section [5] We then present our new 
multi-dimensional Hierarchical Structure Finder (HSF) in 
Section [3l In Section |4l we use our algorithm to detect and 
analyse phase-space structures of a large halo taken from the 
Millennium Simulation. We investigate the space of param- 
eters on which our HSF algorithm depends and try to find 
the best choice of the parameters according to the applica- 
tion under consideration. We also introduce the simulation 
merger tree to follow the evolution of structures in phase- 
space. This allows us to analyse in detail a few representa- 
tive cases. This is followed by a quantitative analysis of HSF 
structures in the space and time domain. We also discuss the 
bimodal nature of the substructure population, in terms of 
bound structures versus tidal tails and tidal streams. Finally, 
in Section[5]we give a summary and present our conclusions. 



2 STRUCTURE FINDERS 

An important step in the analysis of cosmological A'^-body 
simulations is to search for virialized dark matter haloes. 
These are commonly defined as regions around local den- 
sity maxima enclosed by a certain isodensity contour. The 
exact definition of such a border changes from method 
to method. The simplest and the most popular technique 
for findi ng virialized haloe s is the friends-of-friends (FOF) 
method (|Davis et al.l Il985l ) , which links together particles 
which are separated by less than a fixed length b. Usu- 
ally b is set to 0.2 times the mean inter-particle separation, 
which corresponds to finding haloes with overdensity ap- 
proximatel y equal to 178 time s the mean background den- 
sity pmoan ( Colc fc Lacevl[l996l ) . The mass function of haloes 
identified by the FOF method is in good but not perfect 
agreement with the predictions of the Press-Schechter the- 
ory. However, the meth od tends t o link together structures 
across fine bridges (e.g. iLukic et~a l. 2008) and it is not ca- 
pable of detecting substructures inside the virialized haloes 
themselves. A comp arable method is the spherical overden- 
sity algorithm (SO, iLacev fc Collll994l ) which searches for 
local density peaks and then grows around them spheres out 
to a radius where the enclosed mean density satisfies a pre- 
scribed overdensity criterion. By definition, the SO method 
finds only spherical structures. It does not link structures to- 
gether with artificial bridges as FOF does, but it may count 
mass twice in certain cases. 

However, for current high resolution simulations, one 
needs to find not only isolated haloes but also their inter- 
nal substructures. One of the first methods which made it 
possible to find such structures is t he hierarchical friends-of- 
friends scheme (jKlvpin et al.|[l999l ). in which a set of differ- 
ent linking parameters, 6, is used to identify multiple levels 
of substructures inside haloes. 

To distinguish haloes and their substructures in rich 
environments, each detected structure is then usually tested 
against an additional binding criterion. This dynamical cri- 
terion uses information from velocity space to guarantee that 
each structure not only exists but will survive for a longer 
period of time. 

In the spirit of the SO and FOF met hods, the bound 
density maxima (BDM. iKlypin et "aLlll999l ) and DENMAX 
|Ge lb fc B ert schingej Il994 ) methods were proposed. In 
BDM, particles are grouped in spheres around local den- 
sity maxima and are then progressively unbound. In DEN- 
MAX, particles are grouped together when they converge to 
the same local density maximum if they are moved along 
local gradients, calculated on a rectangular grid. In a fi- 
nal additional step they are attached to groups identified 
by the FOF method and then their total binding energy 
is checked. This method was generalised in the SKID al- 



gorithm (iGovernato et allll997l ). in which the local density 
and its gradient are calculated directly at the particles posi- 
tions with the SPH method. A similar b ut simpler method 
was implemented in the HOP algorithm IjEisenstein fc Hug 
Il998h . in which each particle is connected to the one with 
the highest density (found by SPH) among its A^ngb clos- 
est neighbours (with A'ngb ranging typically between 10 and 
20). In this way, space is divided into peak-patches that are 
then combined into the final structures. 

The HOP method gave rise to new structure finders 
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Figure 1. Sketch of the hierarchical structure finder. SI, S2, S3 and S4 are four different structures found by the HSF algorithm. We 
start to grow structures from local maxima which are marked with (a). We grow them down by looking at local neighbours with higher 
density (b), and we connect particles properly to structures. Particles on the border of two structures are always connected to the larger 
one. When we find a saddle point (c) we first apply our Poisson noise criterion. When the structure is not significant (S4), we connect all 
its particles to the most massive partner (cl). If both structures are significant, we apply to them the cut or grow criterion, and when 
one structure is significantly less massive than its companion, we grow only the most massive one (c2), or if structures have comparable 
masses, we grow both of them (c3). 



such as SUBFIND, ADAPTAHOP, VOBOZ and PSB. The 
differences between them are sometimes quite subtle. In 
SUBFIND (|SpringeI et al.ll200ll '). each particle marks its two 
closest neighbours with higher density among the A'ngb clos- 
est ones. Then particles are sorted by SPH-density in de- 
scending order. Particles without higher density neighbour 
are marked as the centres of new structures (local density 
maxima) . Then the structures grow down in density till they 
reach border particles, called saddle points, which have two 
higher density neighbours belonging to two different struc- 
tures. The smallest structure is marked as a structure candi- 
date and then both structures are joined together. Structure 
candidates are arranged in a hierarchical tree and are suc- 
cessively unbound, going from bottom of the tree to the top. 
Particles which are not bound to a structure are attached 
to its larger parent structure. 

Even thou gh ADAPTAHOP 

jAubert. Pichon fc Colombil |2004 ) constructs the tree 
of structures in the opposite way to SUBFIND, from 
bottom to top, the main ideas are very similar. First, it 
grows peak-patches around local density maxima as in HOP 
and then finds border particles and saddle points among 
them. In addition to SUBFIND, each structure is checked 
against Poisson noise to infer its level of significance. In 
ADAPTAHOP, contrary to SUBFIND, only particles above 
saddl e points define structures. VOBOZ (jNevrinck et al.l 
I2OO5I ) on the other hand uses Delaunay tessellation to define 



particle densities and neighbourhood relations, and also 
checks the significance level of structures agai nst a specific 
Poisson noise criterion. The PSB algorithm (iKim fc ParkI 
I2OO6I ) uses a grid as in DENMAX to find local density 
maxima and saddle points, but then constructs a hierar- 
chical structure tree in the same way as in SUBFIND. In 
PSB, particles below saddle points are first attached to all 
structures which are above them and then are assigned to 
individual structures following the process of unbinding. In 
addition to the standard unbinding procedure, PSB takes 
into account tidal force criteria. 

Even though many of the above structure finders use 
velocity information for the purpose of a gravitational un- 
binding procedure, none of them uses the full six dimen- 
sional phase-space information. However, the advantage of 
such an approach is that structures can be defined in a much 
more natural way in phase-space. In particular, they have a 
higher contrast than in position space. In fact, many struc- 
tures such as streams and caustics are well defined only in 
phase-space. 6D FOE l|Diemand et al.ll2006l l is the first im- 
plementation of a structure finder working directly in phase- 
space. It is conceptually a simple extension of EOF based 
on a six dimensional distance measure, using a fixed global 
scaling between position and velocity space. The proposed 
method finds only local phase-space density maxima and 
then grows spheres around them like in BDM algorithm. 

In this paper, we propose a new universal multi- 
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dimensional Hierarchical Structure Finder (HSF) which is 
used here to find phase-space structures in cosmological N- 
body simulations. The algorithm employs, in a higher num- 
ber of dimensions, a similar approach to SUBFIND and 
ADAPTAHOP, but with a new and very eff'ective cut or 
grow criterion, controlled by a connectivity parameter a, to 
separate accurately structures from each other. 



3 THE HIERARCHICAL STRUCTURE 
FINDER (HSF) 

Our goal is to find the hierarchy of dark matter haloes and 
subhaloes, which are defined by locally overdense regions in 
phase-space. The main difference between our approach and 
previous structure finders is that we focus on all kind of 
phase-space structures even those which are not self-bound, 
such as tidal streams. To enable comparisons, we however 
implement, in addition to our base algorithm, also an un- 
binding step. The HSF can be run on a full simulation to de- 
tect all the haloes and the subhaloes population, or on stan- 
dard groups found by a FOF algorithm with, e.g., b = 0.2. 

Prior to the identification of structures, our HSF al- 
gorithm estimates the local phase-space density and the 
local phase-space neighbourhood of each particle i n the 
sample. Following the proposal of iMacieiewski et al.l (|2008| ) 
for optimum local phase-space density estimation, we use 
the SPH method with A'sph neighbours found in the 
adaptive metric computed by the EnBiD algorithrrQof 
ISharma and Steinmet j l|2005l ). We performed a small mod- 
ification of the EnBiD algorithm to make possible the out- 
put for each particle of the phase-space neighbourhood with 
A'ngb closest neighbours among the A^aph (in the proper lo- 
cal adaptive metric frame). It is worth mentioning that both 
phase-space density and neighbourhood estimations by the 
EnBiD algorithm are computationally inexpensive and are 
almost as fast as standard three-dimensional SPH estima- 
tors. 

We define phase-space structures as the regions grown 
around local density maxima by following the local density 
gradient. To find such structures, HSF uses a modified ver- 
sion of SUBFIND, which redistributes particles below saddle 
points in a new fashion. In the first step, the HSF algorithm 
finds locally overdense regions in phase-space by tracing iso- 
density contours identified by saddle points. In addition, we 
test on each saddle point if structures are statistically sig- 
nificant when compared to Poisson noise as in ADAPTA- 
HOP. Particles below phase-space isodensity contours can 
in principle be attached to many structures simultaneously 
but our aim is to attach each of them to only one structure. 
To do that, we use a simple but robust cut or grow crite- 
rion depending on a connectivity parameter a, which allows 
us to reconstruct a multilevel hierarchy of structures within 
structures. In our implementation, each saddle point defines 
a connecting bridge between two structures. Particles be- 
low this saddle point are attached to one of the structures, 
or are redistributed between both of them according to the 
structure masses. We here consider two different cases: 

• the two structures have comparable masses: we attach 



SPH-AM in the notation of IMacieiewski et al.l ||2008|) 



particles below the isodensity contour to both structures in 
the same manner as above the isocontour. Border particles 
are attached to the most massive structure; 

• one structure is significantly less massive than the other 
one: we mark it as a structure and attach all particles below 
it only to the most massive one. 

While we apply the cut or grow criterion, we create a hierar- 
chical binary tree of structures by connecting each structure 
to its more massive partner if one exists. This way of cutting 
works like a second Poisson noise criterion and it allows one 
to grow only structures which are significant. 

In detail, the HSF algorithm, sketched in Figure [TJ 
works as follows: 

(i) For each particle, we estimate the local phase-space 
density with SPH-AM and the local adaptive metric en- 
vironment using Nngh neighbours. We usually perform the 
SPH interpolation with A^aph ~ 64. This value represents 
a good compromise between filtering of Poisson noise and 
identification of faintest significant structures. We find that 
the final results are rather insensitive to the choice of Angb- 
Our favourite value is Angb = 20, similar to what is used 
with HOP, ADAPTAHOP and SUBFIND. Then for each 
particle, we find the set A of its neighbours among the Angb 
which have higher density than the particle. We sort the set 
A ascendingly according to local neighbourhood distances 
(closest particles are in the beginning of the list). Then we 
take the two closest elements of A and put them in a second 
set B. This set can be empty or contain one or two elements. 

(ii) We sort the particles by decreasing phase-space den- 
sity and, following this ordering, we attach each particle to 
different structures according to the following rules: 

(a) The set B is empty. This means that the parti- 
cle does not have any neighbour with higher density: we 
found a local maximum and we mark the particle as the 
beginning of a new structure. 

(b) The set B contains one or two particles which be- 
long to the same structure: we attach the particle to this 
structure; or set B contains two particles which belong 
to different structures Sm and S„, and the Sm structure 
is marked as a more massive partner of 5'„: particle is 
attached to Sm (this is a border particle). 

(c) The set B contains two particles which belong to 
different structures Sm and S„ and the structures are not 
on each other's more massive partner list: it means that we 
found a saddle point and we perform the marking Sm > 

Sn, Sn ^ 5'm, Or Sm — Sn. 

The way this marking is performed in detail can be 
described as follows: 

(1) First, we check the level of significance of struc- 
tures Sm and Sn when compar ed to Poisson noise 
(|Aubert, Pichon fc Colonibill2004l ). Let (Sm) and (Sn) 
be the first and the second structure's average density 
and Psaddie be the density of the saddle point connecting 
them. Each structure is significant if 



{S)> psaddlc 



1 



(1) 



where P is the "/Ja" level of significance of the structure 
(in our tests /3 is set between and 4), and N is the 
number of particles belonging to the structure. If one 
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of the structures is not significant, then we attach all of 
its particles to the second structure. In the case where 
both structures are not significant, we attach all the 
particles to the structure which has the highest maxi- 
mum density. 

(2) If both structures are significant compared to 
Poisson noise, we test them against the cut or grow 
criterion. Let \Sm\ and \S„\ be the masses of our struc- 
tures up to this saddle point and \Sm\ > \Sn\, then we 
mark structure Sm as more massive partner of Sn- If 
ISmla > with a €]0, 1], then structure Sn is more 
than 1/q times less massive than Sm and we attach all 
the particles below this saddle point to Sm- 



SUBFIND unbound 



(3) If -jg^ G]a, i[, we consider that both structures 
have the same order of mass: we attach the saddle point 
to the most massive structure and all particles below are 
attached according to the rules we set before. 

(iii) Finally, a structure containing less than A^cut parti- 
cles is considered insignificant, and all its particles are at- 
tached to its more massive partner. If a structure with less 
than A'^cut particles does not have a more massive partner, 
we put it on the list of fuzzy particles. 

(iv) At the end of this process, we obtain a hierarchical 
tree of structures. Each particle belongs only to one struc- 
ture or to the background (fuzzy list). In addition, we add 
to our algorithm a final step in which we check each struc- 
ture against an unbinding criterion. Once we have marked 
its more massive partner for each structure, we sort them re- 
cursively such that the larger partners (parents) are always 
after the smaller ones (children). Then we unbind structure 
after structure from children to parents and add unbound 
particles to the larger partner. For each individual structure, 
we calculate the gravitational potential. We set the struc- 
ture centre as the position of the particle with the minimum 
potential and the velocity centre as the mean velocity. We 
calculate the kinetic energy of each particle relative to the 
mean velocity of the structure. All the particles with posi- 
tive total energy are marked and, in that ensemble, 1/4 of 
the ones with positive total energy are removed. We repeat 
this process iteratively (starting with a new gravitational 
potential calculation) up to the moment when we stay with 
bound particles only. If the structure has less than A'cut par- 
ticles after the unbinding process, then we mark it as not 
bound and attach all its particles to its more massive part- 
ner or put them on the fuzzy particles list. To speed up the 
calculation of the gravitational potent ial, we use th e tree 
algorithm implemented in GADGET-2 (|Springelll2005h . 

Most halo finders such as DENMAX, BDM, SKID, 
SUBFIND, ADAPTAHOP and VOBOZ use a two step pro- 
cedure for finding the structures. First, they assign as many 
particles as possible to each individual structure in three 
dimensional space by tracing local overdensities (Figure [21 
top left panel). When we move to phase-space diagram (Fig- 
ure [2] top right panel), we however immediately observe 
that there are many particles belonging to different velocity 
structures. The unbinding process (Figured second row of 
panels) then cleans up all these spurious velocity structures. 
In the four bottom panels, one can observe the results ob- 
tained with the six dimensional HSF algorithm. This method 
allows us to attach particles to structures in a more natu- 
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Figure 2. Appearance of our Millennium simulation halo (colour 
ranging from green to red, scaling logarithmically with phase- 
space density) and superposed to it, one of its largest substruc- 
tures found by different algorithms (grey pattern). Left panels: x-y 
position space; right panels: radius r - radial velocity Vr phase- 
space. Prom top to bottom, the grey pattern corresponds to the 
substructure found respectively by (i) SUBFIND before unbind- 
ing, (ii) SUBFIND after unbinding, (iii) HSF before unbinding, 
(iv) HSF after unbinding. 



ral way, because it treats both position and velocity space 
(Figure [2l third row of panels). Note that, after unbinding, 
the structures detected by the HSF algorithm are more ex- 
tended than with standard algorithms working in position 
space (Figure [21 bottom panels), an indication that more of 
the mass belonging to the substructures is recovered. 
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Figure 3. Mass distribution of the substructures as a function of 
the ratio between the substructure mass and the mass of the main 
halo. Here, the influence of the choice of the shot noise control 
parameter /3 (Eq. [TJ on the mass profile is tested. 



Figure 4. Mass distribution of the substructures as a function 
of the ratio between the substructure mass and the mass of the 
main halo. Here, the influence of the connectivity parameter a 
and the importance of the unbinding process are tested. 



4 RESULTS FOR A TEST HALO FROM THE 
MILLENNIUM SIMULATION 

To test our algorithm we use a large halo extracte d from 
the Millennium Simulation ijSpringel et al.l l2005al ). The 
main cosmological parameters of this ACDM-simulation are: 
nm = 0.25, h = 0.73, D.A = 0.75 and erg = 0.9 {Ho = 
lOOh km s~^ Mpc~^). The simulation volume is a periodic 
box of size 500 fe~^Mpc and individual particles have a mass 
8.6 X 10* /i~^Mq. In our analysis we take the second largest 
FOF halo at redshift z = 0, which has 3.83 x 10® particles. 

This section is organised as follows. In Section 14.11 we 
discuss the influence of the main parameters in our algo- 
rithm on the results. In Section [4.21 we use the merger tree 
history to follow both qualitatively and quantitatively the 
evolution of structures. In particular, the structures identi- 
fied by HSF are cross-correlated with their counterpart prior 
to merging with the main halo. Finally Section [4.31 studies 
the properties of the substructure population obtained with 
HSF and its bimodality in terms of bound structures versus 
unbound tidal tails and tidal streams. 



4.1 Choice of the main parameters in the 
algorithm 

In the following, we check the influence of the different pa- 
rameters on the structures found by our HSF algorithm. A 
basic parameter setup is given by A^sph — 64, A'^„gb — 20, 
/? = 0, Q = 0.2, iVcut = 20. We adopt the notation 
(-^sph, A'^ngb, /3, Q, dimension, (B)ound/(UN)bound) to label 
each set of parameters. When the dimension is set to 3D, 
we mean the three dimensional position space, whereas 6D 
means six dimensional phase-space. Unless mentioned oth- 
erwise, we use the HSF algorithm without additional un- 



binding step. SUBFIND is in fact one of the versions of our 
algorithm, characterised by the following parameter setup 
(64, 20, 0, always cut smaller partner, 3D, B). In our analy- 
sis, we shall call all the particles in a FOF group a halo, 
the largest structure of the FOF group a main halo, and all 
other structures substructures. 

Figures|3]and|4]present the number of substructures per 
logarithmic mass bin scaled to the main halo mass. While 
testing the different parameter setups, we find that a and 
/3 are the most influential. Figure |3] shows that parameter 
/3 is important for the smallest structures: visual inspection 
suggests that a higher /3 helps to preserve small scale con- 
nectivity, e.g. between tidal tails and the bound component 
of the substructures. The connectivity parameter a has a 
similar effect, but on the whole mass range instead of small 
structures only, as illustrated by Figure [4] Because of the 
partial degeneracy between a and /3, we prefer to use (3 = 
in most of our analyses. The choice of a indeed influences 
connectivity as follows: when a = 0.2, the main halo always 
wins the cut or grow criterion and all structures are cut by 
it; when a = 0.02, the largest substructures can grow in- 
side the main halo; when a = 0.01, all small substructures 
grow more aggressively and the halo is divided into more 
small parts. In brief, focusing on bound structures calls for 
a value of a of the order of 0.2, while if one is interested in 
all substructures including tidal streams, it is better to set 
Q ~ 0.01 — 0.02. In the later case, tuning up (3 can help to 
control the small scale connectivity. 

Using our base parameter setup, we now compare HSF 
bound structures with those given by SUBFIND. The HSF 
algorithm works in six dimensional phase-space, where the 
distribution of particles shows much more contrast than in 
position space alone. Because of that, HSF can better trace 
contours of individual substructures and attach more parti- 
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Figure 5. Ratio of mass of each HSF bound structure divided by 
the mass of its SUBFIND counterpart as a function of SUBFIND 
structure mass. The central curve corresponds to the median value 
of the ratio calculated over 10 logarithmic bins along the x-axis, 
taking into account only bins containing 2 points or more. The 
two additional green curves on each side show 1 cr errors estimated 
from the variance of points in each bin. 



cles to them. Even after the unbinding step, HSF therefore 
attaches more particles to the substructures than SUBFIND. 
This is illustrated by Figure O where the ratio between 
the mass of HSF bound structures and the mass of their 
SUBFIND counterparts is plotted: HSF attaches on average 
~ 1.1 more mass to small structures than SUBFIND and up 
to twice more to the largest ones. 

The left panels of Figure [6] compare bound structures 
found by both methods in position space. The area of each 
circle is proportional to the structure mass. 

With the parameters set up chosen here, HSF finds 
around 10% more structures, mostly small ones, in the out- 
skirts of the main halo and clearly attaches more mass 
than SUBFIND to most of them. Nevertheless, the spa- 
tial distributions of HSF and SUBFIND substructures are 
nearly the same, as expected. To complete this visual in- 
spection, the right panels of Figure [S] compare the 200 
largest bound structures found by both methods. Most 
of the HSF structures are matched by SUBFIND, except 
that they are more extended. As mentioned before, many 
of these structures have 1.1 — 2 times more mass in HSF 
than in SUBFIND. Interestingly, this confirms the mass ex- 
cess found around SUBFIND substructures in a compari- 
son of simulation with gravitational len sing observations by 
iNataraian. De Lucia fc Springel (|2007l '). 

If bound substructures are counted as a function of 
maximum circular velocity instead of mass, a much closer 
agreement is found, however. This is seen in Figure [T] where 
the cumulative velocity functions of bound substructures for 
HSF and SUBFIND are compared. HSF tends to find a few 
more small substructures, but both algorithms essentially 



identify the same set of more massive structures, confirming 
the results above. 



4.2 Phase-space structures and merger tree 

In the following sections we describe a method to follow back 
in time the structures detected by our HSF, paying partic- 
ular attention to the definition of what we use for initial 
halo. Then we study in detail a set of specific but repre- 
sentative cases. The goal of this analysis is to physically 
understand the nature of the structures found by our algo- 
rithm, before and after unbinding. In particular, we aim to 
separate clearly tidal streams from compact bounded sub- 
haloes. With the additional time information, we can also 
associate tidal streams to objects in the stage they were 
prior to merging with the main halo. We can also study 
quantitatively how in general phase-space structures evolve 
in time. 



4-2.1 Evolution with time and the merger tree 

To study in detail the nature of phase-space structures found 
by the HSF algorithm we use the merger tree historjjf] to fol- 
low their evolution backwards in time. Then we count how 
many particles are shared between each structure prior to 
merging with the main halo and its counterpart detected by 
HSF at z — 0. This process uses pieces of information which 
are alr eady available for the processed Millennium Simu- 
lation (jSpringel et al ] l2005bl ) and it is divided into three 
steps: (i) crosscorrelating the HSF structure catalogue with 
the SUBFIND one, (ii) following the evolution of SUBFIND 
structures using the already implemented merger tree his- 
tory and (iii) using each particle's universal indejifl to follow 
structures at different output times. We now explain each of 
these steps in turn. 

(i) cross-correlation between HSF and SUBFIND: Infor- 
mation about structures in the Millennium Simulation is 
organised in terms of two levels: first, particles are attached 
to different FOF groups (found with b = 0.2). Then, in each 
FOF group, they are separated into the main halo, the sub- 
structures found by SUBFIND, and unbound 'fuzzy' parti- 
cles if present. Running the HSF algorithm with the base 
parameter setup (64, 20, 0, 0.2, 6D UB/B) provides a phase- 
space structure list. Then, for each member of that list, the 
SUBFIND substructure sharing the largest possible number 
of particles with it is identified. If the HSF structure shares 
less than 20 particles with any SUBFIND substructure, it 
is put into an unmatched list. In the opposite case, we call 
this SUBFIND substructure a seed of the HSF structure. 



^ Each branch of this tree corresponds to the evolution in time of 
a dark matter halo as a stand alone structure, while each node of 
it corresponds to the event of merging between 2 halos or more. 
Note that, due to the coUisionless nature of dark matter, the halos 
can pass through each other and separate again: in practice, the 
structure of such a tree can be non-trivial. 

^ The universal index of a particle is just a number associated to 
each single particle in order to identify it unambiguously, which 
is useful for analyses of Lagrangian nature such as performed in 
this work. 
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Figure 6. Left side panels: Spatial distribution of HSF bound structures and SUBFIND structures, with at least 20 particles. The area 
of each circle is proportional to the structure mass. Right side panels: First 200 most massive bound substructures identified by HSF and 
their SUBFIND counterparts. Particles belonging to the same substructure share the same colour. 



(ii) Following SUBFIND structures back in time: Infor- 
mation about the time evolution of structures is stored in 
th e Millennium Simula tion in a merger tree (more details 
m ISpringel et al]|2005bD . We use the tree information which 
gives for each halo or substructure its most massive progen- 
itor, if there is any. Once the list of seed SUBFIND sub- 
structures is obtained, each of them is traced back in time 
by following its most massive progenitor recursively up to 
the moment when this past structure was the main halo 
of a FOF group. This is the last occurrence of the struc- 



ture as being distinguishable as an isolated halo. We store 
the redshift of this event and all particles belonging to the 
main halo found in this way are denoted as the initial halo. 
There is a small number of substructures which do not have 
a proper progenitor, they are dropped from the analysis. 

(iii) Using each particle's universal index to follow struc- 
tures at different times: In the last part of the procedure, 
we link together the information gathered in the previous 
two steps. For each HSF structure identified at the present 
time, we find its initial halo and, with the help of universal 
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Figure 7. Cumulative count of substructures as a function of 
maximum circular velocity. We here compare results of HSF for 
bound structures with substructures identified by SUBFIND in 
the same halo. 

indices, we identify shared particles, i.e. particles belonging 
both to the initial and final structures. We carry out ex- 
actly the same analysis for HSF bound structures and for 
SUBFIND itself. 

In addition, during this process, we gather for each sub- 
structure information about the position of its centre and its 
velocity at various times (we use the SUBFIND definition 
for the structure centre), and similarly for the position and 
velocity of the centre of the main halo. With these pieces 
of information at hand, we can define an orbital count by 
determining the number of times a substructure's radial ve- 
locity with respect to the centre of the halo changes sign, 
which each time signals completion of what we call an orbit. 
This definition requires that there are enough snapshots to 
catch orbit details. This is the case for most substructures, 
probably all, although this statement is not easy to check 
rigorously. 

4-2.2 Definition of the initial halo 

In our analysis of the time evolution of the structures, we 
adopt the main haloes of the FOF groups found by SUB- 
FIND as 'initial haloes'. Another possibility is to chose for 
each initial halo all the particles belonging to the FOF 
group. In the first case, the analysis is simplified by the fact 
that we look only for the evolution of one isolated compo- 
nent. However this component represents only one part of 
the halo. In the second case, a FOF group can sometimes 
have a few main halo candidates joined by small artificial 
bridges (up to 20% of FOF groups show such a feature, e.g. 
iKim fc Par^l2006l ) but can in fact be tidally disrupted such 
that its components get away from each other. 

To demonstrate the effects described above, we choose 



one particular structure in our test halo. The top panel of 
Figure [8] shows all the particles belonging to the initial halo 
traced to redshift z = 0, while the second row of panels cor- 
responds to the full traced initial FOF group. This structure 
goes around the main halo one time (its orbit is shown in 
the third row of panels of Figure [5}. The initial FOF group 
is tidally disrupted during this process and its various com- 
ponents are clearly separated from each other. SUBFIND 
recognises the central part of the bound object (bottom pan- 
els of Figure [8]). The HSF bound structure contains more 
particles (fourth row of panels in Figure ^ . These parti- 
cles belong to tidal tails, but are in fact still gravitationally 
linked to the structure. The HSF structure (prior to unbind- 
ing) contains 55% of the particles of the initial halo, and 
reproduces perfectly its shape (on third row of Figure [S| . 

4-2.3 Qualitative analysis of structure evolution 

To better understand all the processes at play during struc- 
ture evolution, we study in greater detail eight different cases 
displayed in Figure [S] The colours in the figure are coded as 
follows: 

• Green and red particles belong to one structure found 
by the HSF algorithm (without unbinding): green particles 
belong to the initial halo, while the red particles do not 
belong to it. 

• Black particles belong to the SUBFIND seed of the HSF 
structure. 

• Blue particles belong to the initial halo, but do not 
belong to the HSF structure. 

• Yellow particles belong to the initial halo and belong to 
any HSF structure, besides the one we take for the current 
analysis. We mark particles in yellow only for structures in 
the last three rows of panels of Figure [9] 

The particles are plotted in the following order: blue, red, 
green, black, and finally yellow. Various structures parame- 
ters are listed in Table [1] for each of the 8 cases considered 
here. The pink curve shows the orbit of the object inside the 
main halo. We now discuss in detail each of these 8 cases. 

(i) The first case corresponds to the largest structure 
found by HSF. It entered the main halo recently, at redshift 
z = 0.17. HSF identifies a large fraction of 84% of the ini- 
tial structure, and unbinding worsens that number by only 
about 3%. On the other hand, SUBFIND identifies only 37% 
of the content of this rather massive structure. 

(ii) The second structure is at a moment just before cross- 
ing for the first time the main halo centre and starts to be 
significantly tidally disrupted. The HSF structure still con- 
tains 97% of particles of the initial main halo, while the 
SUBFIND bound structure accounts for only 26%. Indeed, 
HSF manages to attach to the structure unbound particles 
which already crossed the main halo centre and contribute 
to a tidal tail. 

(iii) The third structure entered the halo at redshift z = 
0.56 and made an orbit inside it. This is the reason why we 
identify only 55% of the initial structure, but still more than 
SUBFIND (35%). HSF attaches some additional particles to 
the structure, i.e. particles that do not belong to the initial 
main halo, but in fact we found that a large fraction of them 
belong to the initial FOF group. 
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Figure 8. Follow up of a particular substructure (in colour) of our Millennium test halo (in grey). From left to right: x—y position 
space, radius r-radial velocity Vr, radius r-phase-space density /. In the two left columns of panels, the colour traces the logarithm of 
the phase-space density (from dark grey to light grey or from dark blue to red). In the right column of panels, the colour just traces 
the projected particle density. Prom top to bottom: Initial main halo traced to z = 0, initial FOF group traced to z = 0, HSF unbound 
structure, HSF bound structure, and SUBFIND structure. 
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Figure 9. Properties of some chosen structures. From left to right: x-y position space diagram, radius i — radial velocity Vr diagram, 
radial velocity tangential velocity vt diagram. Further description in the text (section l4.2.3l l. 



(iv) The fourth case corresponds to a structure that just 
passed the main halo centre. 88% of the initial structure is 
still identified. HSF is capable of connecting to the struc- 
ture (unbound) particles which were strongly affected by 
tidal stripping but still belong to the structure. In this case 



for example, it joins some particles having high negative ve- 
locity. 

(v) The fifth case corresponds to the example of a struc- 
ture that has gone around the main halo centre. HSF identi- 
fies 90% of the initial structure and SUBFIND only 44%. We 
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Figure 9. Continued. 



can clearly see here the typical "z" shape of the structure 
in {r,Vr) space, which correspond to a scaled down version 
of the phase-space diagram of the halo (that we see only in 
half here) . 

(vi) The sixth row of panels corresponds to a rare occur- 
rence when HSF partly fails. The HSF structure contains 



65% of particles of the initial one. It is falling inside the 
main halo centre and some of the particles already crossed 
the centre. Parts of the initial structure are identified as 
other HSF objects: the initial halo shares 72% of its parti- 
cles with all detected HSF substructures. There is however a 
large fraction of particles associated with the HSF structure 
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Table 1. Properties of eight chosen structures. Column description: nr.: structure number, in the same order as listed in the text and 
in Figure O orbits: number of orbits inside the main halo; HSF par.: number of particles belonging to the HSF structure; Initial FOF 
group par.: number of particles belonging to the initial FOF group; Initial halo par.: number of particles belonging to the initial main 
halo; SUBFIND: fraction of the particles of the initial halo that are still identified in SUBFIND (in HSF bound, HSF, any identified 
HSF structures, respectively for the next columns). 



nr. 


orbits 


z 


HSF par. 


Initial FOF group par. 


Initial halo par. 


SUBFIND 


HSF bound 


HSF 


any HSF 


1 





0.17 


78791 


109835 


76567 


37% 


81% 


84% 


94% 


2 





0.09 


37417 


39511 


32383 


26% 


84% 


97% 


99% 


3 


1 


0.56 


36248 


85225 


53186 


35% 


44% 


55% 


57% 


4 


1 


0.24 


35520 


43202 


35463 


26% 


43% 


88% 


91% 


5 


1 


0.32 


33647 


30432 


28604 


44% 


60% 


90% 


92% 


6 


1 


0.62 


17971 


17053 


14364 


32% 


38% 


65% 


72% 


7 


2 


0.83 


14477 


140799 


119621 


2% 


5% 


11% 


12% 


8 


4 


1.50 


854 


10696 


9896 


3% 


6% 


8% 


25% 



(in red) that should just belong to the background. Indeed, 
we checked that most of them cannot even be associated 
with the initial FOF group. 

(vii) The seventh row of panels shows the typical case of 
a massive structure which, after only two orbits (so it passed 
nearby the halo centre only twice), already dissolved in the 
main halo, because of massive tidal disruption. Even though 
the HSF structure contains only 11% of the particles of the 
initial halo, we find that other HSF substructures match 
some parts of the initial halo: 90% of particles inside such 
substructures come from the initial halo, although some of 
them belong to other members of the initial FOF group. 
In other words, it means that our algorithm is capable of 
finding remnants of tidal tails. All blue particles can not be 
distinguished from the main halo anymore and correspond to 
the part of the structure which has been completely diluted. 

(viii) In the last case, we take a structure which merged 
with the main halo at high redshift z — 1.5, and already 
made four orbits inside. As expected, this structure has ex- 
perienced strong tidal striping: the HSF structure contains 
only 8% of the initial structure. Most of it indeed already 
dissolved inside the main halo. Still, some remnants were 
identified as disjoint HSF components and represent 25% of 
the initial structure. 

To conclude this section, objects which recently entered 
the main halo and typically made up to one orbit inside it are 
in most cases fiiUy recovered by HSF. When the structures 
make more orbits, they are more tidally disrupted, especially 
when they come close to the halo centre, so the fraction of 
particles identified decreases. HSF still finds in most cases 
remnants of strongly disrupted objects as individual tidal 
streams detached from the bound component (if this later 
still exists). 

4-.2.4 Quantitative analysis of structure evolution 

To test the performance of HSF quantitatively, one can for 
example study the fraction Msharod /Minitiai of particles in- 
side the initial halo found at the present time in the corre- 
sponding HSF structure, as a function of redshift of merg- 
ing with the main halo (top left panel of Figure llOp or as a 
function of initial mass (top right panel of Figure llip . In- 
deed, one expects a strong correlation between the value of 
Afsharcd/Afinitiai and the initial halo mass and the redshift z 



of merging. Obviously, the higher z, the larger the number 
of orbits (Fig. I12p and the larger the number of particles 
lost due to tidal stripping (top left panel of Figure [10)) . At 
low redshift, z < 0.3, where the number of orbits is typi- 
cally less than one, there is no strong structure evolution 
in phase-space and HSF identifies 80 — 100% of the initial 
structure mass (upper left panel of Fig. I10|l . There are a few 
structures at redshift z ~ 0.1 — 1.0 for which HSF can find 
only a very small fraction of their initial mass. All of them 
are tidal remnants. 

The effect of unbinding on the ratio Maharcd/.Afinitiai is 
shown in the lower left panel of Fig. [11] (as a function of 
mass) and upper right panel of Fig. [TU] (as a function of red- 
shift). Obviously, after unbinding, the fraction of particles 
recovered by HSF decreases significantly, even for a small 
redshift of merging z < 0.3. Indeed, a significant fraction of 
the mass in substructures is contained in tidal tails that are 
very well identified by HSF, at least for z < 0.3, but that 
are not bound any more to the substructures. Note that, 
as expected, the SUBFIND (bound) substructures are not 
very different from the HSF bound ones, except that they 
contain a slightly smaller fraction of the mass of the initial 
structures (lower right panel of Fig. Illjl . 

Another important test of our structure finder consists 
in examining the fraction of particles inside each HSF struc- 
ture that is shared with the initial halo as a function of 
initial halo mass (top left panel of Figure llip . HSF finds 
that for massive objects, 80 — 90% of the present structure 
mass belongs to the initial halo. For smaller structures, the 
scatter is higher and only around 40% of the particles found 
in HSF objects belong to initial halos. This actually means 
that for many small structures, the HSF algorithm attaches 
more particles than they had before. This it mainly due to 
the fact that we consider for the initial stage only the main 
part of the halo and not its substructures: disentangling sub- 
structures from the main halo remains an ambiguous pro- 
cess, and structures identified at the present time can con- 
tain part of the mass of the substructures inside the initial 
halo. In the process, we also do not take into account parti- 
cles surrounding the initial halo which were not selected by 
the FOF (yet) but were still infalling onto our Millennium 
test halo and participate to its background density. As a re- 
sult, additional particles can be associated to the final HSF 
structure, and this effect is expected to become stronger if 
the mass of the identified structure is small. 
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Figure 10. Mass loss of structures as a function of tlie redsliift z of merging with the main halo. Top-left panel: the mass loss is 
presented as the ratio MsjiarodZ-'^^initiali where Afgijarod is the mass in common between the HSF structure and its counterpart (of mass 
^initial) just prior to merging with the main halo. Top-right panel: same as top left one, but for bound HSF structures. Bottom-left 
and bottom-right panels: same as the top ones, but the mass loss is presented as the ratio between total final mass and initial mass, for 
HSF bound structures and SUBFIND (bound) structures, respectively. On all the panels, the symbol size is proportional to Mjnitiai- In 
addition we plot the median value (in red) and the a errors calculated in 10 logarithmic bins (in green), with at least 10 structures per 
bin. 



Keeping that in mind, we can now study the ratio 
Mx/Afinitiai between the total mass of the identified bound 
structure and the mass of the initial halo as a function of 
redshift of merging, independently of whether the particles 
are shared between initial and final stage or not, as shown 
in the lower panels of Fig.[TO]for HSF (X=HSF bound) and 
SUBFIND (X^SUBFIND). The resuhs are not very differ- 
ent for both codes, except that, as already extensively ar- 
gued previously, this ratio is slightly larger for HSF than for 



SUBFIND. And because of the effect just discussed above 
(top left panel of Fig. Ilip . Mx/Minitiai can be larger than 
unity at low redshift. 

While the correlation between the mass loss due to tidal 
stripping and the redshift of merging is quite obvious, the 
relationship between mass loss and initial mass is less evi- 
dent, given the limited amount of statistical occurrences we 
analyse here (one single large cluster-sized Millennium halo). 
Figure[TT] (right panel and bottom panels) indicates that the 
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Figure 11. The nia,ss, A^sj^a^j-Q{j , in commoii between structures detected at present time by HSF and SUBPIND and their counterpart 
-of mass Minitial— just prior to merging with the main halo, is studied in a fractional way as a function of Afinitiai- Top-left: ratio between 
Afgijared ^'nd the mass of the HSF unbound structure. Top-right, bottom-left and bottom-right panel: ratio between Mgijared ^^id Mi^itiai, 
respectively for HSF unbound, HSF bound and SUBFIND structures. The symbol size is proportional to redshift z of merging with the 
main halo. In addition we plot the median values (in red) and the cr errors, calculated in 10 logarithmic bins (in green), with at least 10 
structures per bin. 



mass loss is significantly smaller for initially light structures 
than for initially massive ones. To demonstrate that unam- 
biguously, we perform a more accurate analysis of the global 
mass loss displayed on the lower panels of Fig. 1101 by fitting 
analytically the redshift and the mass dependence. To do 
so, we divide the initial mass of the structures into four log- 
arithmic bins. Then for each bin, we fit the mass loss as a 
function of redshift in logarithmic coordinates (Figure I13p 
with the following convenient parametric form 



Jl^ = I . (2) 

Minitial iz/Zs)'^il-{- Z/Zsp 

The best fitting parameters, found by a standard least- 
square method, are listed in Tabled The fact that mass loss 
is more pronounced for more massive objects is clear, and 
was also to be expected. This behaviour can simply be ex- 
plained as follows: small structures are more st rongly bound, 
because they are more concentrated (e.g. lAngulo et al.l 
200^), so they do not lose as much mass as large structures 
from tidal stripping. Indeed, the largest structures are less 
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Figure 12. Number of orbits each structure made inside the main halo as a function of redshift of merging of the structure with the 
halo. The symbol size is proportional to M^^n^^i, expressed in units of MQh~^ . In addition, we plot the median value (in red) and the 
a errors (in green) calculated in 10 logarithmic bins with at least 10 structures in each bin (for convenience, binning is performed on y 
axis). In our sample the structures do not made more than 6 orbits inside the main halo, before they disappear. 
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compact and are more sensitive to dynamical friction. As a 
result, they are strongly disrupted while they are orbiting 
around the main halo. They also tend to have more radial 
orbits. 



4.3 Bound subhaloes, tidal tails and tidal streams 

By studying the merger tree history, we could show that the 
HSF algorithm is capable of finding both substructures and 
their tidal tails. As we noticed, some of the tidal tails are 
still connected to their host substructures, while others are 
recognised as separate objects. We now study this bimodal- 
ity more carefully. 

To better separate tidal streams from the bound coun- 
terpart of substructures, we now take a small value of the 
connectivity parameter, a — 0.01. In the following analysis, 
we shall study the five following populations: 

(i) Bound substructures found by HSF; 

(ii) Bound substructures found by SUBFIND; 

(iii) All HSF substructures (before unbinding); 

(iv) Unbound HSF structures: substructures found by 
HSF which disappear during the unbinding process, such 
as tidal streams; 

(v) Bound HSF structures: substructures found by HSF 
(along with their tidal tails) which remain after the unbind- 
ing process. 

To analyse substructures properties, we estimate, for each of 
them, the phase-space density maximum /max and the min- 
imum /min. These 2 quantities are measured for the full set 
of particles belonging to the substructure prior to the un- 
binding process. To better measure high phase-space density 



peaks, we use the SPH-AM EnBiD method with 32 neigh- 
bours. 

In the top-right panel of Figure 1141 the distribution of 
the/max values of is shown. It is clearly bimodal, and this 
property is in fact independent of a. It is straightforward 
to understand the origin of the bimodality. The local phase- 
space density can for in stance be approximated as follows 
(|Macieiewski et aLllioOsI ): 

f(x V) - c::p[ ["-""W]' ] (3) 

Two cases can the be considered. In the centre of a bound 
substructure, i.e. a standalone structure that survived the 
unbinding process, the local density p is high and the lo- 
cal velocity dispersion (t(x) is low, which gives a high local 
phase-space density maximum. On the other hand, when 
substructures are disrupted by strong tidal forces, their lo- 
cal density p decreases and their velocity dispersion, cr(x), 
increases, so their peak phase-space density is lower. 

These statements can be directly checked by unbinding 
the substructures found by HSF. On the top-right panel of 
Figure 1141 the high density maxima peak of the distribu- 
tion is dominated by the bound substructures, as expected. 
There is a small fraction of bound substructures for which 
/max resides in the lower density maxima regime. We checked 
that this happens only for the smallest substructures with 
around 20 particles, for which Poisson noise fluctuations 
start to be significant. The lower peak of the distribution 
of values of /max is mainly occupied by unbound substruc- 
tures. There are still some unbound substructures residing 
in the higher peak. They have less than 100 particles and can 
be considered as "Poisson clusters" (even in Poisson noise it 
is possible to find high density contrasts). 

Note that the high /max distribution peak is very sharp. 
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Figure 13. Mass loss as a function of redshift of merging z for masses binned in 4 logarithmic bins (dotted curves) with its smooth fit 
given by equation ^ (thick curves). For each dotted curve, the number of bins is equal to 2\/N, where A'^ is the number of samples. 
These two panels are equivalent to bottom panels of Fig. 1101 To make adequate fitting, we perform Levenberg-Marquardt least-square 
minimisation with sigma errors set from Poisson noise counting distribution. Masses are expressed in units of Mq/i"^. 



Table 2. Parameters used in Eq. jJl to fit results presented in Figure [T3l 



mass min 


mass max 


HSF Xs 


HSF rj 


HSF 7 


SUBFIND 


SUBFIND 77 


SUBFIND 7 


1.7 X 10"' 


9.8 X 10"' 


0.54 


0.04 


0.78 


0.30 


-0.00 


0.66 


9.9 X lO^" 


5.6 X 10" 


1.02 


0.05 


1.82 


0.40 


0.00 


1.06 


5.6 X 10" 


3.1 X 10^2 


20.11 


0.00 


28.31 


1.46 


-0.03 


3.01 


3.2 X 10^2 


1.8 X 10" 


19188.89 


0.01 


28851.01 


18.93 


-0.07 


25.61 
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Figure 14. Disentangling tidal streams from bound substructures: phase-space density distributions. Top-left: distribution function of 
logarithm of phase-space density / estimated for all particles belonging to each category of substructures as indicated inside the panel 
(100 logarithmic bins). The black dashed lines represent the best fitted Gaussian functions for the main halo found by SUBFIND, the 
main halo found by HSF and the unbound structures found by HSF (including the tidal tails of bound structures) . This means in fact 
that for each of these components, / is lognormal if the fit is good. To make adequate fitting, we perform Levenberg-Marquardt least- 
square minimisation with sigma errors set from Poisson Noise counting distribution. Top-right, bottom-right: distribution function of 
substructures maxima, /max, and minima, /min (50 logarithmic bins). The substructures are separated into unbound components (blue) 
and bound ones (green), while the black curve corresponds to all the substructures. Bottom-left: distribution function of substructures 
phase-space "peakness", defined as Cf = /max//min- 



Table 3. For each substructure we measure its maximum phase-space density /max, minimal value /mim maximum 3D density pmax, 
minimum one Pmim ^^•^ phase-space density "peakness", = /maxZ/mim and 3D density "peakness", Cp = Pmax /Pmin- Phase-space 
density is quoted in M0/i^kpc~'^km~''s^, while p is expressed in units of total average density (p). 



Structure class 


/max 


/mill 


Cf 


Pm 


iX 


pmin 


Cp 


HSF unbound 


4.0 X 10^5 


3.3 X 


10-6 


8.7 


2.1 X 


10^ 


69.2 


6.1 


HSF bound 


0.2 


2.7 X 


10-*^ 


1.1 X 10-'^ 


1.4 X 


10* 


31.2 


35.5 
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Figure 15. Disentangling tidal streams from bound substructures: projected 3D density distributions. Top-left: distribution function 
of normalised density 1 + 5 = p/{p) estimated for all particles belonging to each category of substructures as indicated on the panel 
(100 logarithmic bins). Top-right, bottom-right: distribution function of substructure /(p), minima Pmin/ (p) (50 logarithmic 

bins); bottom-left: distribution function of substructures density "peakness" defined as Cp = Pmax/Pmin- 



correspon ding to /max — 0.2M(^/t^k pc~'^km~'^s'^. As already 
noticed in lMacieiewski et al.l (l2008h . all the bound substruc- 
tures pre sent approximately the same value of /max (see also 
IVass et al. 2008). This property could be simply an upper 
bound imposed by numerical resolution or set by the dy- 
namics, or more likely a combination of both (e.g., iBinnevI 
l2004l : IVass et al.ll2008l ). The second peak, dominated by tidal 
streams, is less pronounced, although still quite well defined, 
with a maximum at /max — 4.0 x 10~^MQ^^kpc~^km~^s^, 
a value about 3.7 orders of magnitude lower than what is 
found for bound structures (all the values are summarised 
in table [31). This shows again the very clear separation be- 
tween bound structures and tidal streams. 

Another way of separating various substructure popula- 



tions consists of measuring their "peakness", i.e. the param- 
eter Cf — /max//min, whete /min is the minimum value of 
the phase-space distribution function of the HSF structures 
(prior to unbinding) . The advantage of the peakness param- 
eter is that, as opposed to /max, it does not depend on a 
specific choice of units: a structure with a bound component 
should present a peakness parameter very large compared to 
unity, contrary to a pure tidal stream. The measurement of 
Cf is however meaningful only if /min is well defined. This is a 
priori not obvious as one expects /min to be very sensitive to 
local fluctuations in the noise, which indeed affect the local 
topology strongly. We checked that in fact /min is a robust 
statistics, as suggested by the rather symmetric behaviour 
of the curves shown in the bottom-right panel of Figure 1141 
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Table 4. Mass in each substructure class compared to the total 
mass in our Millennium test halo. 



Structure class mass 

SUBFIND bound 12.4% 

HSF bound 18.8% 

HSF unbound a = 0.2 22.5% 

HSF unbound a = 0.01 31.6% 

HSF a = 0.2 41.2% 

HSF a = 0.01 50.4% 



The distribution of measured values of c/ is shown in the 
bottom-left panel, and presents of course the same bimodal 
nature as /max- For instance, one finds that c/ is typically 
of the order of 10^ for bound structures, while it is only of 
the order of 10 for tidal streams. 

Finally, the top-left panel shows the distribution of mea- 
sured values of / for each dark matter particle. In this plot, 
particles left over after unbinding HSF substructures, i.e. 
belonging to the tidal tails of these substructures, are put 
on the list of unbound substructures. The high phase-space 
density region is dominated by bound substructures, which 
is consistent with the observations we made for the /max dis- 
tribution function. Note that HSF bound substructures are 
more extended into lower phase-space density regions than 
SUBFIND ones and are more likely to overlap in terms of 
density with unbound streams. There is in total almost 19% 
of mass in HSF bound substructures to compare with 12.4% 
in SUBFIND ones (see Tabled}. This additional mass in 
HSF bound substructures comes from particles which were 
not found with the saddle point algorithm working in 3D po- 
sition space on which SUBFIND is based. This means that 
the total bound mass of substructures strongly depends on 
the cutting criterion applied to the 3D density field p, even 
with the additional unbinding procedure. 

An examination of the top left panel of Figure [TJ] 
suggests that it is possible to perform an optimal cut 
on /, with a value chosen between 3 x 10~^ and 3 x 
10"'* MQ/i'^kpc~'^km~''s'^ so that most particles with phase- 
space density above this threshold belong to bound sub- 
structures. Such a criteri on was used b efore in the literature 
to mark substructures l|Stadel et al.l '2008) . Tidal streams 
and possibly signatures of caustics occupy the middle range 
of phase-space densities, with 31.6% (22.5% for a = 0.2) 
of the total FOF halo mass belonging to them, which is 
more than for bound substructures. Similarly as for bound 
substructures, we can set some lower limit around 10~^ on 
the phase-space density and claim that most particles with 
higher value of / than this limit belong to substructures of 
some kind (bound or unbound) . The low phase-space density 
regime is indeed dominated by the main halo component. 

As a final note on the upper- left panel of Figure 1141 we 
found that the shape of the distribution function of values 
/ observed for each component has interesting properties: 
it is very well fit by a lognormal distribution both for the 
main halo component found by SUBFIND and HSF, and the 
unbound subst ructures found b y HSF. This complements 
the findings of IVass et al.l (|2008t ). who performed a similar 
analysis but used a more add-hoc approach to separate vari- 
ous components contributing to the phase-space distribution 



Table 5. Best parameters of the Gaussians fitted to the distribu- 
tion function of the logarithm of phase-space density estimated 
for all particles belonging to each category of substructures (top- 
left panel of Figure [T4t . 



Structure class 


mean 




cr 


error 


SUBFIND main halo 


4.7 X 


10- 


-6 


0.73 


7.39 


HSF main halo 


2.7 X 


10- 


-6 


0.55 


1.11 


HSF unbound substructures 


1.4 X 


10- 


-5 


0.65 


0.34 



function. The best fitting parameters of a Gaussian on the 
logarithm of / are given in Table [S] The interpretation of 
these results did not seam straightforward to us, so we de- 
cided to leave it for future work. Certainly, a path to follow 
is to examine the arguments developed by I Coles fc JonesI 
(|l99lf ) to explain the close to lognormal behaviour of the 
projected 3D density, p, relying on the continuity equation 
and the positivity of the density. 

In practice, in observations of, e.g.. X-rays clusters or 
gravitational lensing, the 3D density p (or its projection on 
the sky) is usually used to model dark matter haloes in- 
stead of the phase-space density /. To illustrate how the 
previous results translate in terms of p. Figure [15] is similar 
to Figure [141 but the calculations are performed for the nor- 
malised density 1 + 5 = p/{p) instead of /. The 3D density 
is measured using EnBiD's SPH kernel with 32 neighbours. 
Contrary to Figure 1141 the distribution function of values 
of Pmax (upper right panel) shows only one peak. The dif- 
ference between bound and unbound structures shows much 
less contrast (see Table [3] for numerical estimates of typi- 
cal values of pmax, Pmin and Cp = pmax/pmin). lu particular, 
bound structures present a large spread on their 3D local 
density maxima of about 2 orders of magnitudes, in contrast 
with what happens with /max, and they are more difficult to 
disentangle from their unbound counterpart, even with the 
peakness, Cp, although this latter quantity seems to have a 
better separating power than pmax (lower left panel). 

Interestingly, the particle density distribution diagram 
(top-left panel of Figure I15p is populated in a different 
way from what happens in phase-space. In particular, tidal 
streams occupy the low density region although they still 
spread over a large dynamic range, while the main halo 
dominates the high density regime. Bound substructures 
are rather subdominant and spread over the whole dynamic 
range. 

To complete this section. Figure [T6l and Figure [TTl show 
the appearance of bound structures, unbound ones, and of 
the smooth part of the halo after removal of all HSF struc- 
tures. There is a subtle but significant difference between the 
two figures. In Figure [16] the top panels show only the bound 
part of the bound substructures, while the top panels of Fig- 
ure[T7]show the bound structures along with their tidal tails. 
In the middle panels of Figure [16] particles both belonging to 
unbound structures and particles removed from the bound 
structures during the unbinding process are shown. In con- 
trast, the bottom panels of Figure [TTl show only particles 
belonging to unbound structures. This results in an asym- 
metry in middle right panel of Figure 1161 which reflects the 
fact that structures passed through (or nearby) the centre of 
the halo one more time in the upper part of the phase-space 
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Figure 16. Appearance of bound and unbound structures in our Millennium test halo. Top panels: particles belonging to HSF bound 
structures (so unbound particles are removed). Middle panels: particles which belong to HSF unbound substructures or are left over in 
the tails of HSF bound substructures after the unbinding process. Bottom panels: the main halo after removal of all substructures. From 
left to right; x—y position space, radius i — radial velocity Vr phase-space. The pictures are computed in 3 steps as follows: (i) division of 
space into a three-dimensional equally spaced grid with N = 400 divisions across each x,y, z axes, (ii) calculation of the mean density / 
of all particles inside each cell and (iii) projection of this density on the x — y plane by taking in each z column the cell with the highest 
density. To enhance the contrasts, equalisation of the histograms in log / was implemented. 
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Figure 17. Same as Figure 1161 but tlie set up is sliglitly diilerent: in the top panels, the HSF structures which are seeds of bound 
structures are displayed entirely, including their tidal tails. In the bottom panels, only HSF structures which are unbound are displayed. 



diagram than in the lower part. Tidal disruption is indeed 
more significant and thus removes particles with higher value 
of / in the upper part of the phase-space diagram than in 
the lower part. Not surprisingly, the asymmetry disappears 
in the bottom right panel of Figure 1171 Note that bound 
structures are absent in the region close to the main halo 
centre, as expected. Note as well the rather elongated tidal 
streams, in particular close to the halo centre, in the bottom 
right panel of Figure [T7] and the middle right panel of Fig- 
ure 1161 These are the left overs of structures disrupted by 
strong tidal forces. In the bottom right panel of Figure 1161 
the main halo still presents, after cleaning, some filamentary 
structures, which are parts of tidal tails, or less likely, sig- 
natures of caustics. It can be cleaned even more by using a 
smaller value of the connectivity parameter a. 

5 DISCUSSIONS AND CONCLUSIONS 

We introduced a new universal multi-dimensional hierarchi- 
cal structure finder (HSF) which was employed here to study 



dark matter structures in six dimensional phase-space. The 
algorithm used, for each particle, the phase-space density 
and local neighbourhood estimated with the SPH method 
with adaptive metric implem ented in the EnBiD package 
(jSharma and Steinmetd |2005| ). To detect structures, HSF 
builds on the SUBFIND and ADAPTAHOP algorithms, 
with the introduction of a new simple but robust cut or 
grow criterion depending on a single connectivity parameter 
a. 

The main steps of the algorithm are as follows: (i) lo- 
cal phase-space density maxima are detected and structures 
around them are grown by following local gradients up to 
saddle points; (ii) at each saddle point level, the density / of 
each structure is compared to Poisson noise: the structure is 
kept if / is /3 times more significant than the Poisson r.m.s. 
noise level. At the same time its mass and the mass of its 
partner (connected to it through the saddle point) are mea- 
sured and the cut or grow criterion is applied: if one struc- 
ture is a times smaller than its neighbour, then all particles 
below the saddle point are attached to the neighbour. When 
the two structures have comparable mass within a factor a. 
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they are both set to grow down as before. This criterion al- 
lows us to better trace substructures in phase-space, with 
a good control of the effect of Poisson noise, which is very 
important in this rather sparsely sampled space. 

We demonstrated the potential of HSF on a large FOF 
dark matter halo taken from the Millennium Simulation. 
Our tests show that /3 and especially a are important con- 
trol parameters. To better study the smallest possible struc- 
tures, /3 should be set close to 0. The smaller a, the more 
subtle the structures found by the algorithm. In our analysis, 
we give preference to a = 0.2, which provides a good bal- 
ance between finding the finest possible substructures and 
not overgrowing them. This value of a is particularly ap- 
propriate when an additional binding step is performed. In 
contrast, an analysis of tidal tails is best carried out with 
small a, around 0.01 — 0.001, which separates structures into 
smaller pieces. It is possible to use it in combination with 
/3 = 4 — 10, which tends to reconnect the structures together 
in a consistent way, to reconstruct tidal tails rather well. A 
more advanced method of reconnecting phase-space struc- 
tures, by using the topology of the hierarchical tree created 
by the HSF algorithm is under investigation. 

We used the M illennium Simulation merger tree 
ijSpringel et al.1 l2005bl ) to compare the HSF phase-space 
structures found at the present time with the same struc- 
tures traced back to the time just before they enter the 
main halo. While the best three dimensional algorithms used 
presently, such as SUBFIND, manage to find only the main 
part of bound structures, HSF is capable of finding more 
extended bound components along with their tidal tails. 
There is much more information about structure evolution 
still stored in phase-space than in 3D and this information 
can be potentially fully recovered from the data by a six 
dimensional algorithm such as HSF. 

The main results of our analysis in time and space do- 
main are the following: 

• HSF structures contain on average 80 — 100% of the 
mass inside the initial structures up to a redshift of merging 
z = 0.3 — 0.4. This value drops down to 50% for z = 1. On 
the other hand, bound HSF structures contain on average 
80 — 100% of the mass inside the initial halos only up to 
z = 0.09 and 50% up to z = 0.6. This shift in the mass loss 
is caused by the existence of tidal tails, which are joined 
to HSF structures, but do not belong to their bound part. 
In other words we can say that HSF is able to reconstruct 
in most cases the full dynamical structures which enter the 
halo at redshifts as high as z = 0.3 — 0.4. 

• The distribution function of the phase-space den- 
sity maxima /max of HSF structures clearly shows a bi- 
modality. We can explain it by partitioning the struc- 
tures into two distinct groups. In the first group, cor- 
responding to the high phase-space density peak regime, 
/max ~ O.2M0/i'^kpc~^km~'^s^, with a small spread around 
that value, there are mostly bound structures. In the 
second group, corresponding to a 3 orders of magni- 
tude smaller phase-space density regime, /max ~ 3.3 x 
10~^MQ/i'^kpc~'^km~^s'', and a larger spread around this 
value, there are all unbound structures i.e. tidal tails, 
streams and possibly some caustics. In terms of "peakness" , 
Cf = /max//min, whcrc /min is minimum value of / in each 
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Figure 18. Cumulative mass in substructures found by SUB- 
FIND, bound HSF method, and HSF with different choices of the 
connectivity parameter a. 

substructure, this translates into Cf — 1.1 x 10^,8.7 for 
bound and unbound structures, respectively. 

• We noticed, similarly as _Vass_et_al. (2008), that the dis- 
tribution function of the values of / around each dark matter 
particle is close to lognormal for the smooth component of 
the halo and the unbound part of the substructures (tidal 
streams) . 

• We found that there is more mass in bound HSF sub- 
structures than in SUBFIND ones. Figure [TS] shows the cu- 
mulative mass of substructures divided by the total halo 
mass as a function of substructure mass. Around 18.5% of 
halo mass is stored in bound HSF structures, and this quan- 
tity almost does not depend on a. In comparison, about 
12.4% of the mass is stored in SUBFIND bound structures. 
The additional mass in HSF bound structures comes mainly 
from the fact that subhaloes are better defined in phase- 
space and are more extended. However, the set of identified 
bound substructures is nearly identical in both methods, 
and hence the cumulative abundance of substructures as a 
function of maximum circular velocity is the same as well. 

We note that in our test halo, 41.2% of the mass belongs 
to substructures for a = 0.2, 50.4% for a = 0.01, and 55.2% 
for a — 0.001. When we substract from these numbers the 
contribution of bound structures, we find that 22.9% — 36 A% 
of halo mass is stored in unbound structures. This should 
be taken into account when analytical models of halos with 
substructures are proposed. 

• While we would need a larger statistical sample of halos 
to perform robust measurements, we noticed that at fixed 
redshift of merging with the main halo, small structures tend 
to lose less mass than larger ones, in agreement with expec- 
tations based on the higher concentration of smaller halos. 
Also, we found a strong correlation between mass loss and 
the number of orbits a substructure can make inside the 
main halo. 
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When we observe our own Galaxy, we do not have ac- 
cess to different "snapshots" anymore, in stark difference 
with the world of simulations. Instead, we have to be content 
with the data at the present time. However, because we now 
know that our phase-space structure finder can identify dy- 
namical structures that were bound before tidal disruption, 
it can provide totally new insights about the past dynamical 
history of our Galaxy. Within the hierarchical framework, we 
expect that our Galaxy should be made through the merging 
of more than about 100 smaller subcomponents. Comparing 
structures in observational data and simulations can be one 
of the best tests for the theory of hierarchical galaxy for- 
mation, and provide important constraints on cosmological 
models such as ACDM. 

Up to now, we studied only the evolution of dark mat- 
ter, but we can also similarly study the evolution of baryons 
in gas and stars. Galaxies are observed in many different 
ways, ranging from star distributions, velocity and chemical 
properties, to HI measurements, etc. Our phase-space struc- 
ture finder with local metric fitting is in fact implemented 
in such a way that it can be used in any number of dimen- 
sions, where each dimension can have completely different 
physical units. So it is in principle straightforward to use it 
for studying galaxy structure evolution in multi-dimensional 
space with the appropriate probabilistic weightings to take 
into account the noise and holes (missing measurements) in 
the data hypercube. We think that such an approach can 
yield a deeper understanding of galaxy evolution, and looks 
especially promising in light of the upcoming GAIA mission 
l|Gilmore et al.lll998l ) which plans to map the positions of 
around one billion stars in our Galaxy. 
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